Majorana Modes in Driven-Dissipative Atomic Superfluids With Zero Chern Number 
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We investigate dissipation-induced p-wave paired states of fermions in two dimensions and show 
the existence of spatially separated Majorana zero modes in a phase with vanishing Chern number. 
We construct an explicit and natural model of a dissipative vortex that traps a single of these modes, 
and establish its topological origin by mapping the problem to a chiral one-dimensional wire where 
we observe a non-equilibrium topological phase transition characterized by an abrupt change of a 
topological invariant (winding number). We show that the existence of a single Majorana zero mode 
in the vortex core is intimately tied to the dissipative nature of our model. Engineered dissipation 
opens up possibilities for experimentally realizing such states with no Hamiltonian counterpart. 



The search for topological phases of matter in which el- 
ementary excitations exhibit non-Abelian statistics have 
brought two-dimensional (2D) p-wave paired superfluids 
and superconductors to the forefront of theoretical and 
experimental condensed-matter research [1-5]. The bulk 
of these systems is fundamentally intriguing in that it 
reveals physics beyond the Landau paradigm: different 
phases are characterized by distinct values of a nonlocal, 
topological order parameter known as the Chern num- 
ber, and phase transitions occur whenever the topology 
changes, signaled by discontinuities in this integer-valued 
topological invariant. In topologically non-trivial phases 
corresponding to an odd Chern number, vortices with 
odd vorticity have been predicted to carry unpaired Ma- 
jorana fermions, and to exhibit non-Abelian exchange 
statistics as a result [6]. 

In this Letter, we explore the concept of topologi- 
cal order and its connection to the edge physics in a 
non-equilibrium scenario based on engineered dissipa- 
tion [7, 8]. Prior work on quantum-state engineering 
in driven-dissipative systems has shown that topologi- 
cally non-trivial states of many-body Hamiltonians can 
also be prepared as steady states of a dissipative dynam- 
ics [9]. In contrast, we demonstrate that dissipation can 
lead to a novel manifestation of topological order with 
no Hamiltonian counterpart. Specifically, we show that 
spatially separated Majorana zero modes (MZMs) can 
be obtained in a 2D, dissipation-induced p-wave paired 
phase of spin-polarized fermions with vanishing Chern 
number. Remarkably, a phase whose topological nature 
is seemingly trivial -according to the standard diagnostic 
tool provided by the Chern number- can therefore exhibit 
phenomcnological features characteristic of a non-trivial 
one, which ultimately leads to the counterintuitive fact 
that vortices with odd vorticity may obey non-Abelian 
exchange statistics in a bulk with zero Chern number. 

We demonstrate these results in a simple model moti- 
vated by an implementation scheme based on cold atoms 
and optical vortex imprinting [10], where fermion parity 
is microscopically conserved. We show that they hold 



over an extended parameter range, in which we identify 
a topologically non-trivial phase missed by Chern num- 
ber considerations. Critical points are revealed upon in- 
troducing a vortex, in which case we establish the phe- 
nomenology of a non-equilibrium topological phase tran- 
sition characterized by (i) a discontinuity in a topologi- 
cal invariant (winding number) , (ii) divergent length and 
time scales [8, 11-13], and (iii) a divergent localization 
length associated with a MZM bound to the vortex core. 

The basic mechanism behind our findings relies on 
the fact that the introduction of a (dissipative) vortex 
changes the system in two crucial respects: it modifies its 
topology, as argued in Refs. [14, 15], and imposes specific 
(dissipative) boundary conditions. Here we show that our 
model of a vortex with odd vorticity can be mapped to 
a ID chiral fermion problem [16] characterized by a non- 
trivial topological invariant (winding number) = 2 
despite a vanishing bulk Chern number. In such a situa- 
tion, bulk-edge correspondence arguments [4, 17-19] sug- 
gest the existence of a pair of MZMs in the vortex core. 
However, owing to the dissipative boundary conditions 
imposed by the geometry of the vortex core alone -which 
underpins the universal nature of our findings- a single 
MZM only is found in the core. This phenomenon cru- 
cially relies on dissipation, and therefore has no Hamilto- 
nian counterpart. It shows that the potentially harmful 
effect of dissipation on MZMs [20-22] need not be en- 
tirely destructive, but may instead give rise to intriguing 
novel effects. 

Dissipative framework - We consider a system of N 
fermionic sites a\ , ai evolving under a purely dissipative 
dynamics governed by a Lindblad master equation 

N 

d tP = K j2{ L *p L l-H L \Li,p}), a) 

i=l 

where p is the system density matrix, k the damping 
rate, and Li are Lindblad operators which are linear in 
the fermionic operators. The steady state of such dynam- 
ics is pure if and only if the Lindblad operators form a 
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set of anticommuting operators, i.e. {Li,Lj} — for all 
i,j = l,...,N. If so, it can be identified with the ground 
state of the parent Hamiltonian H paront = J2i-^l^i- Al- 
though purity is not required for the existence of topolog- 
ical order [9], we will only encounter steady states that 
are pure in the bulk, whose (bulk) topological proper- 
ties can be inferred from H paient alone. In our dissipa- 
tive setting, the analog of a gap in a Hamiltonian spec- 
trum is a dissipative gap in the Liouvillian spectrum, 
which dynamically isolates the subspaces corresponding 
to the bulk and edge modes, thereby providing the coun- 
terpart of a gap protection through the quantum Zeno 
effect [23] . Most importantly, the counterpart of topolog- 
ical ground-state degeneracy is the existence of a nonlocal 
decoherence-free subspace associated with zero-damping 
Majorana modes 7 = 7^ which satisfy the orthogonal- 
ity condition {^,7} = = for all i (see Sup- 
plementary Information (SI)). Clearly, this condition is 
more restrictive than the one for a zero-energy Majo- 
rana mode of a Hamiltonian, which reads (for -ff par ent) 

[^parent, 7] = E^H^} ~ {4^} L i) = °- ^ Cru- 

cial difference stems from the first "recycling" term on 
the right-hand side of Eq. (1), and is the reason why dis- 
sipation can crucially modify the Majorana physics found 
in the Hamiltonian context. 

The model - We consider a square-lattice system 
driven by so-called "cross" Lindblad operators defined as 
the following quasi-local linear superposition of fermionic 
creation and annihilation operators a\,ai\ 

U = C\+ aj+Ai = pa\ + (al + a\ 2 + a{ } + a\j 
+ ae l4 '(ai 1 + ia i2 - a i3 - ia i4 ), 

(2) 

where f3 £ R, a > 0, <fi € [0, 2ir), and i\, 12, «3 and 14 are 
the four clockwise-ordered nearest-neighboring sites of i. 
The creation and annihilation parts of Li have s- and p- 
wave symmetries, respectively. The dissipative dynamics 
generated by such Lindblad operators can be obtained, 
for example, as the long-time limit of a microscopically 
number-conserving (quartic) dissipative dynamics gener- 
ating phase-locked paired states (see SI). In that con- 
text, the global relative phase <p between the creation 
and annihilation parts of Lj emerges through sponta- 
neous breaking of the global U(l) symmetry, and the 
relative strength a is determined by the average particle 
number. The dimensionless parameter /3, on the other 
hand, can be used to tune the system across phase tran- 
sitions, as will be shown below. 

The steady-state bulk properties of the system arc 
most easily revealed in the infinite-size limit. Defining 
coefficients Uij, Vij such that Lj = ^2j(uijOj +VijOj), the 
momentum-space cross Lindblad operators take the form 
Lk = itkflk + Wka^ k , where Uk, ^k are the Fourier trans- 
forms of , Vij . Properly normalized, they become Bo- 
goliubov quasiparticle operators associated with damping 



rates «k = K A/"k, where A/k = with = (ukjfk)- 
One can easily verify that the cross Lindblad operators 
form a complete set of anticommuting operators, so that 
the system is driven into a unique and pure complex p- 
wave paired state |fi) defined by the condition L^\D,) = 
(for all k) and fully characterized by the real vector 
nk = £ k c£k, where er is a vector of Pauli matrices. Since 

(i) this steady state is pure (i.e. |nt| = I for all k) and 

(ii) time-reversal symmetry is broken due to the complex 
p-wave nature of the state [24], the topological invariant 
relevant to that case coincides with the Chern number 
i^2D commonly used in 2D Hamiltonian systems (see e.g. 
Ref. [25]). Here we find that the Chern number van- 
ishes [26], namely, 

^2D = 5f / d 2 k n k • (ds^nk x <9 fey n k ) = (3) 

JBZ 

(where BZ stands for "Brillouin zone") for all values of 
/? except at the isolated points j3 = and j3 — ±4 where 
the dissipative gap closes (see SI). Since there is no ex- 
tended parameter range with non-trivial topological or- 
der, one naively expects topological features such as iso- 
lated MZMs to be absent when edges or vortices are in- 
troduced. We will show below that such conclusions are 
premature: single, unpaired MZMs are generally found 
in the parameter range < \/3\ < 4 when dissipative 
vortices with odd vorticity are introduced. 

Physically, the special values j3 = 0, ±4 appear as crit- 
ical points since the dissipative gap closes at these values, 
leading to divergent length and time scales characteris- 
tic of a second-order phase transition [12, 13]. However, 
the symmetry of the steady state (encoded in nk) and the 
value of the topological invariant i/ 2 d both are identical in 
the neighborhood of those points. It thus seems as if the 
apparent critical behavior can neither be traced to a con- 
ventional phase transition (with broken symmetry), nor 
to a topological one (with a discontinuity in the topolog- 
ical invariant). Below we will show that the introduction 
of a vortex makes it possible to define another topolog- 
ical invariant. This will allow us to identify (3 = ±4 as 
genuine critical points associated with topological phase 
transitions. 

Introducing a dissipative vortex - We now introduce 
a dissipative vortex, by modifying the annihilation part 
Ai = J2j u ij a j °f tnc above cross Lindblad operators. 
More specifically, we replace the translation-invariant 
coefficients u^ = Uij considered so far by position- 
dependent coefficients of the form Uij = f(rj)e~ llV: >Uij, 
where (rj , <pj ) are polar coordinates defined with respect 
to a particular site io chosen as the center of the vortex 
core. Two crucial ingredients appear in this definition: 
(i) a real, rotation-invariant vortex profile function f(r) 
which goes to zero in the vortex core, and (ii) a vortex 
phase that winds i times around the origin i [l being the 
vorticity). These properties of a dissipative vortex nat- 
urally arise in an implementation with ultracold atoms 
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based on optical vortex imprinting (see SI). They are 
fully analogous to those of a vortex in a Hamiltonian 
scenario. In particular, a 7r-fhix is optically imprinted 
onto the matter system in the case t = 1. 

In what follows, we show that an odd dissipative vortex 
(i.e. with odd vorticity) generically traps a single, iso- 
lated MZM despite the seemingly trivial topological na- 
ture of the bulk. We focus on the simplest case of a single 
vortex with I = 1 and proceed in two main steps. First, 
we demonstrate that the vortex phase crucially modifies 
the topology of the system and construct a mapping to an 
effective ID model. Second, we examine the effect of the 
dissipative "boundary" conditions imposed in the vortex 
core by the profile function f(r), and show that the latter 
are responsible for the existence of a single MZM in the 
core. 

Mapping to a chiral ID wire - To exemplify how an odd 
dissipative vortex changes the topology of the system, 
we consider the case £ = 1 and assume, without loss of 
generality, that the vortex profile function f(r) vanishes 
as r — > and satisfies f(r) = 1 for r > r c . We identify the 
region r < r c as the vortex core. In order to capture the 
generic properties of the bulk, we first focus on an annular 
region centered around the vortex core (Fig. 1(a), left) 
where f(r) is constant and the relevant p-wave operator 
carries the vortex phase, e~ ltp (d x +id y ) (in the continuum 
limit, for simplicity). As made explicit in the SI, our 
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FIG. 1. (Color online) . Mechanism ensuring the existence of a 
single MZM in the core of an odd dissipative vortex, (a), Left: 
1=1 vortex on the plane, characterized by a core (blue) and 
a phase factor premultiplying the p-wave operator d x + id y . 
Right: Mapping from the (grey) annular region around the 
vortex core to the cylinder -where the vortex phase disappears 
from the relevant p-wave operator- and reduction to a (chiral) 
ID wire problem (red) . (b) Topological phase diagram for the 
planar model of Eq. (2) with no vortex (blue, characterized by 
the Chern number V2d) and with a single £ = 1 vortex (red, 
characterized by the winding number i>m)- (c) Dissipative 
boundary conditions of the ID wire inherited from the vortex 
core profile (see text). 



model is, in this region, formally equivalent to a cylinder 
model with p-wave operator of the form d x > +id y i ((x' , y') 
being Cartesian coordinates on the cylinder; see Fig. 1(a), 
right) as in the original translation-invariant model on 
the plane (see Eq. (2)). There exists therefore a one-to- 
one correspondence between the original planar model 
with an £ = 1 vortex and the same model in cylinder 
geometry with no vortex. Physically, this stems from 
the fact that the gauge field e~ tv imposed on the plane 
to describe the vortex naturally arises on the cylinder 
owing to the extrinsic curvature of the latter. 

In order to further extract the essence of the single- 
vortex problem, we revert to the original lattice de- 
scription and take advantage of the translation invari- 
ance along the ^/'-direction. Defining Fourier-transformed 
Lindblad operators L x >(k y >) oc ^ , e y l^x' ,y' i we re- 
duce the system to a stack of ID wires that can be inves- 
tigated along the lines of Ref . [9] . In the two momentum 
sectors corresponding to k y i — and tt, the relevant ID 
Lindblad operators take the form 

U = /3' a\ + (aj_ x + af +1 ) + (-Oi_i + a i+1 ), (4) 

where i indexes the lattice sites in the ^'-direction (such 
that L, = L x >(k v > = or tt)) and /3' = /3 + 2 (j3 - 2) 
for k y > — (w). In these two particular sectors, the sys- 
tem thus reduces to a ID wire with chiral symmetry [18]. 
The steady-state bulk properties of this wire can be un- 
veiled in the infinite-size limit, in which case the Lindblad 
operators form a complete set of anticommuting opera- 
tors, leading to a pure steady state described by a real 
unit vector n^. as above (see Eq. (3)). Owing to chiral 
symmetry, this state can be characterized by a "winding 
number" topological invariant z^ld [27]. As detailed in 
the SI, we obtain 

"id = / dk a ' ( n fe x d k n k ) = 2 (5) 
Jbz 

for |/3' | < 2 (i.e. for < |/3| < 4), and Vxd = otherwise 
(a being a unit vector orthogonal to whose existence 
is guaranteed by chiral symmetry). We thus find non- 
trivial topological order in the parameter range delim- 
ited by the special points /3 = and ±4 where the Chern 
number 1/2D exhibits discontinuities. Fig. 1(b) illustrates 
the resulting topological phase diagram: non-equilibrium 
topological quantum phase transitions occur at the criti- 
cal values /3 = ±4, while (3 = corresponds to an isolated 
point separating two topologically equivalent phases. In 
agreement with the full 2D model, the dissipative gap of 
the ID wire closes at each of these values (see SI). 

Dissipative boundary conditions - The fact that we 
can identify a non-trivial topological invariant in an ef- 
fective ID model describing the region surrounding the 
vortex core strongly supports the existence of interest- 
ing "edge" physics inside the latter. The bulk steady 
state of the ID wire, which is pure, can equivalcntly 
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FIG. 2. (Color online). Numerical results for a single dissipa- 
tive vortex with I = 1 on a square lattice of 35 x 35 sites with 
unit spacing and k = 1. Left: Typical form of the MZMs lo- 
calized in the vortex core and on the edge, respectively (here 
for (9 = 3), and localization length scale £ associated with 
the MZM trapped in the core. Right: Low-lying part of the 
damping and purity spectra for ft — 2, both featuring a gap 
with two zero eigenvalues. All results were obtained for a 
vortex profile f(r) = (r/d)e" (r/ « )2 with d = 10 and £ = 20. 



be described as the ground state of the parent Hamil- 
tonian flparent = Bulk-edge correspondence 

arguments based on H paicnt would therefore suggest the 
existence of z^ld = 2 MZMs at the edges of the wire -i.e., 
in particular, in the vortex core. In our general dissi- 
pative setting, however, bulk-edge correspondence argu- 
ments can only be formulated in the presence of a dissi- 
pative gap and in the absence of modes corresponding to 
subspaces in which the steady state is completely mixed, 
which we refer to as purity zero modes (see SI) . Keeping 
in mind that the presence of dissipative boundary condi- 
tions can potentially give rise to such modes, we now pro- 
ceed to examine the "edge" physics that emerges in the 
vortex core as a result of the dissipative boundary con- 
ditions imposed by the vortex profile function f(r). To 
this end, we first extend the above mapping by reducing 
the inner radius of the annular region shown in Fig. 1(a) 
to zero. The resulting extended ID wire is depicted in 
Fig. 1(c); its first, leftmost site corresponds to the center 
of the vortex core where f(r) vanishes, as shown in blue. 
The fact that f(r) varies from site to site in the vortex 
core crucially leads to a violation of the purity condition 
{Li,Lj} — (for all i,j) which is satisfied in the bulk, 
as mentioned above. One can easily verify that the anti- 
commutator {Li,Lj} increasingly deviates from zero for 
\i — j\ < 2 upon approaching the left edge of the wire. As 
a consequence, the steady state increasingly loses purity 
and departs from the ground state of -ff pa rent featuring 
a pair of MZMs. Remarkably, one and only one MZM 
survives at the edge of the wire -or, equivalently, in the 
vortex core- in the full parameter range < \f3\ < 4 as- 
sociated with a non-trivial topological phase. This mode 



is explicitly constructed in the SI, and is shown to be 
exponentially localized, on a characteristic length scale 
£ = ~ l/llog (1^1/2)1 which diverges at |/3'| = 2, 

i.e. at (9 = ±4. At these values which coincide with the 
critical points found above, the norm of the wavefunction 
associated with the MZM diverges with the length of the 
wire, while the dissipative gap closes in the bulk. This 
further confirms the onset of a non-equilibrium topologi- 
cal quantum phase transition. The second MZM naively 
expected in the vortex core, by contrast, acquires a finite 
damping rate and becomes an environmental degree of 
freedom with no correlations with the rest of the system. 
As shown in the SI, any zero mode of -Hparent which ac- 
quires a finite damping rate is effectively traced out of the 
system in steady state, independently of the initial con- 
ditions. We refer to such a mode with no Hamiltonian 
counterpart as an intrinsic purity zero mode, in accor- 
dance with the fact that its existence is intrinsic to the 
dissipative dynamics itself. 

The above findings are supported by extensive numeri- 
cal simulations, confirming the existence of a single MZM 
trapped in the vortex core in the full parameter range 
< |/3 1 < 4 for arbitrary vortex profile functions, as well 
as the divergence of the corresponding localization length 
upon approaching the critical points (3 = ±4 (see Fig. 2). 
As expected, an intrinsic purity zero mode is found in 
the vortex core. Similar features are obtained for odd 
vortices with I > 1; even vortices, in contrast, do not ex- 
hibit any MZM, as expected from the fact that the vortex 
phase can be gauged away in that case (see e.g. Ref. [1]). 
Our numerical results generally support the main conclu- 
sion that odd dissipative vortices generically trap single, 
isolated MZMs despite the seemingly trivial topological 
nature of the system. Following the arguments of Ref. [9] , 
we expect such vortices to exhibit non-Abelian exchange 
statistics when braided around each other through adia- 
batic parameter changes. 
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DISSIPATIVE FRAMEWORK 

In the following section, we discuss the properties of 
the dissipative dynamics considered in the main text and 
its relation to (dissipative) Majorana modes. We start by 
briefly reviewing a mathematical framework originally in- 
troduced in Refs. [1, 2] which provides an alternative to 
the second-quantized formalism used in the main text. In 
this description -made possible by the quadratic nature 
of the dynamics- the density matrix is replaced by a cor- 
relation matrix T which encodes the single-particle cor- 
relations and fully describes the system when Gaussian 
initial states are assumed. The information contained in 
the Lindblad operators is encoded in a damping matrix 
X and a fluctuation matrix Y. As we will show, the ma- 
trix Y shares the formal properties of a Hamiltonian in 
a first-quantized (matrix) representation [3] . The steady 
state, however, is determined by the combined properties 
of X and Y. Dissipative Majorana zero modes, in partic- 
ular, can be identified as zero modes of X. We will give 
precise conditions for the existence of such modes. 

A crucial difference between the ground state of a 
Hamiltonian and the steady state of a dissipative dynam- 
ics lies in the purity of the state. A proper assessment 
of purity can be made by examining the purity spectrum, 
i.e. the spectrum of the matrix L 2 ; pure states are sig- 
naled by a flat (degenerate) purity spectrum with eigen- 
values equal to — 1, while completely mixed subspaces 
are indicated by zero eigenvalues which can be traced to 
the existence of purity zero modes. In the Hamiltonian 
ground-state setting, the state is pure by construction, 
up to a possible degenerate subspace whose purity is not 
determined by the dynamics itself. Such a degenerate 
subspace is found for example in the presence of Majo- 
rana zero modes. In that case, starting from an initial 
high-temperature state, the purity of the Majorana sub- 
space generally vanishes. Since the mixedness is extrinsic 
to the dynamics, and rather results from specific initial 
conditions, we refer to the Majorana zero modes as ex- 
trinsic purity zero modes. The dissipative analog of this 
situation is given by dissipative Majorana zero modes. 
Using the quantum optics terminology, such modes can 
be identified as giving rise to a decoherence-free subspace, 
i.e. to a subspace that is decoupled from dissipation. In 
the dissipative setting, however, a different situation can 
be encountered: the purity of the Majorana subspace can 



be pinned to zero by virtue of the dissipative dynamics, 
independently of the initial conditions. In that case, we 
refer to the corresponding Majorana zero modes as in- 
trinsic purity zero modes. Below we establish criteria for 
the occurrence of this phenomenon with no Hamiltonian 
counterpart. We also derive useful necessary and suffi- 
cient conditions for pure steady states. 

Majorana zero modes and associated 
decoherence-free subspace 

We recall the form of the dissipative dynamics consid- 
ered in the main text, namely, 

n 

d tP = KY J (UpL\-\{L\L i ,p}), (1) 

i=l 

with Lindblad operators Li that are linear combinations 
of fermionic operators a^,aj (j = 1, . . . , N). As shown in 
Refs. [1, 2], the fact that Eq. (1) is quadratic allows for 
a convenient simplification: introducing a basis of Majo- 
rana operators C2j_i = al+a i} c 2 i = \{a\—ai) and assum- 
ing Gaussian initial states, one can focus exclusively on 
the time evolution of the real antisymmetric correlation 
m,atrix Tij = ^ tr([cj, Cj]p) which fully describes the state 
of the system. Expressing the Lindblad operators in the 
vector form Lj = if c, where c T = (ci, c 2 , . . . , c 2 n), one 
can recast Eq. (1) as the following matrix equation: 

d t T = -{X,T} + Y, (2) 

where X = 2 Re M and Y = 4 Im M, with M = J2i h ® 4 
(note that we have set k = 1 in Eq. (1); this factor can be 
absorbed in the definition of Lj, without loss of general- 
ity). Physically, the real symmetric matrix X describes 
the damping; the real antisymmetric matrix Y, on the 
other hand, describes the fluctuations in steady state T 
through the steady-state equation {X, T} = Y. We thus 
refer to X and Y as the damping and fluctuation matri- 
ces, respectively. Their explicit form reads 

Y =-2* EIU (1^4-1^ If)- 

The damping matrix is real, symmetric, and positive 
semidefinite (by construction). It can therefore be spec- 
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trally decomposed as X — Y^j=i ® V J) with eigen- 

values Xj £ ]R> and associated eigenvectors Vj € R 2Ar . 
The real eigenvectors Vj form an orthonormal basis and 
define a natural basis of Majorana operators jj = vjc = 
7J (7 = 1,2,..., 2iV) associated with X (we remark that 
this Majorana basis is not necessarily local, as opposed 
to the original basis given by c T = (a, c 2 , . . . , c 2 n ) )- The 
fact that they are real implies that l\vj = (if v.,)* for all 
j, such that 

n 

X jk ee vjXv k = 2 Rc {(If v,)(lf v fc )*} = Xj S jk (4) 

i=l 

with Aj- = 2££ =1 |lfv j | 2 , and 

n 

y jfe = vJy Vfc = 4^Im{(lf Vj )(lfv fe )*}. (5) 

1=1 

Physically, the eigenvalues Xj correspond to the damp- 
ing rates associated with the Majorana modes jj. As 
shown above, they are solely determined by the amount 
of overlap if Vj between the eigenvectors Vj of X and the 
vectors L_ representing the Lindblad operators in the orig- 
inal (local) Majorana basis. Since lfvj = {Lj,7j}, they 
can alternatively be seen as determined by the amount 
by which the Majorana modes jj anticommute with the 
Lindblad operators Li. In general, we will refer to a 
Majorana mode -jj as a Majorana zero mode whenever 
the associated damping rate Xj vanishes. In light of the 
above discussion, a unit vector v e R 2N defines a Ma- 
jorana zero mode 7 = v T c if and only if the following 
equivalent conditions are satisfied: 

(i) Iv = 0; 

(ii) If v = for all i e {1,2, ...,n}; 

(iii) {Li, 7} = for all i E {1,2,..., n}. 

Remembering the explicit form of the matrix Y , these 
conditions also imply Yv = (the converse is generally 
not true). 

The existence of m < 2N Majorana zero modes gener- 
ally implies the existence of a (m-mmod 2)-dimcnsional 
subspace Vq that is left invariant under the dissipative dy- 
namics of Eq. (2). This stems from the fact that the state 
space of the system is isomorphic to the space of all real 
2 N x 2N antisymmetric correlation matrices T (assuming 
Gaussian states), and is therefore even-dimensional (with 
dimension 2N(2N — l)/2), as must be Vq. As a conse- 
quence, Majorana modes can only decouple from the dis- 
sipative dynamics in pairs. In the original Hilbert-space 
description of the system, Vq translates as a decoherence- 
free subspace [4] of dimension 2( m_mmod2 ). 



Purity zero modes and robustness of the Majorana 
zero modes 

In the matrix representation defined above, the purity 
of a state with correlation matrix T can be inferred from 
the spectrum of T 2 , which we refer to as the purity spec- 
trum. One can easily convince oneself that T describes 
a pure state if and only if P 2 = —1, i.e. iff the purity 
spectrum is flat with all eigenvalues equal to —1. In par- 
ticular, a mode 7 = w T c corresponding to an eigenvector 
w of F 2 with eigenvalue zero can be identified as a purity 
zero mode, i.e. as a mode associated with a subspace in 
which the state is completely mixed. We distinguish two 
types of such modes: (i) the intrinsic purity zero modes, 
which appear in steady state independently of the ini- 
tial conditions, and (ii) the extrinsic purity zero modes, 
which solely result from mixed initial conditions and thus 
disappear if the system is prepared in a pure initial state. 
(Note, however, that a completely mixed initial state is 
a rather generic assumption when a spatially separated 
pair of Majorana modes is present.) 

Keeping the above considerations in mind, we now pro- 
ceed to examine the effect of dissipative perturbations 
on Majorana zero modes. We consider the general case 
where the dissipative dynamics gives rise torn < 2N Ma- 
jorana zero modes j a = v^c, and define V ' C R 2N as 
the m-dimensional subspace associated with these modes 
(with Greek indices denoting quantities pertaining to the 
latter). In addition, we assume that (Gaussian) corre- 
lations are initially present between the Majorana zero 
modes, as determined by the initial correlation matrix 
F a p(t = 0) = I tr([7 Q , Jp]p(t — 0)). By definition, the 
dissipative dynamics restricted to V ' is trivial, namely 
dt^ap = for all pairs of indices a, j3 € {1,2, ... ,m}. 
In order to examine under what conditions the evolution 
remains trivial in V ' when the system is perturbed, we 
extend the dissipative dynamics by introducing an addi- 
tional Lindblad operator L n+ \ = lf +1 c. The X and Y 
matrices restricted to V ' -which were identically zero- 
then become 

X aP ee vf Av„ = 2 Re{(lf +1 v Q )(lf +1 v^n, (6) 
Y ap ee vlYvp = 4 Im{(lf +1 v a )(lf +1 v^)*}. (7) 

Assuming that L„ +1 acts quasi-locally and that the Ma- 
jorana zero modes are localized and spaced sufficiently 
far away from each other (such that {^ ai 7/3} = v a v P = 
up to exponentially small corrections), two cases can oc- 
cur: cither (i) lf +1 v Q = {L n+ i,j a } = for all a e 
{1,2, ...,m}, or (ii) lf +1 v M = {L n+1 ,^} ^ for a sin- 
gle (i€ {1,2,..., to}. Clearly, the dissipative dynamics 
restricted to V ' remains trivial in the first case: the Ma- 
jorana zero modes survive, and all correlations initially 
present in V Q ' are preserved over arbitrary long times. 
In the second case, however, the situation changes: al- 
though the off-diagonal elements of X and Y remain triv- 
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ial owing to the local nature of L n+ i and "f a -such that 
Y remains identically zero- the diagonal element X^ of 
the symmetric matrix X becomes non-zero. More specifi- 
cally, the Majorana mode 7^ acquires a non-zero damping 
rate A M = 2|lJ +1 v M | 2 (independent of the system size), 
and the dissipativc dynamics reads d t T a ^ = — A^F^ 
for all a G {1,2,..., m}. In other words, the introduc- 
tion of an additional quasi- local Lindblad operator L n+ \ 
which anticommutes with all Majorana zero modes but 
7 M has the effect of destroying all non-local correlations 
initially present in Vq that involve the latter. After a 
time t ~ 1/A M , 7^ is effectively traced out of the sys- 
tem and emerges as an intrinsic purity zero mode: for 
arbitrary initial states, a zero eigenvalue appears in the 
steady-state purity spectrum. 

In summary, we have shown that quasi-local dissipa- 
tive processes can generally remove Majorana zero modes 
from the system, thereby resulting in a loss of purity. A 
necessary and sufficient condition for a specific Majorana 
zero mode 7 to disappear under the effect of a particu- 
lar dissipative process given by a Lindblad operator L is 
given by {L, 7} 7^ 0. Although we have used the word 
"perturbation" above, we remark that this condition can 
be satisfied in a controlled way, by carefully engineering 
dissipation. In fact, the model presented in the main 
text provides us with a paradigmatic example of such 
controlled dissipation: due to the generic form of the en- 
gineered dissipative boundary conditions, the removal of 
a single Majorana zero mode naturally occurs (see main 
text). 

Hamiltonian ground states vs. Liouvillian steady 
states 

The parent Hamiltonian H parcnt = XL ^l-^i intro- 
duced in the main text plays a crucial role in understand- 
ing the similarities and differences between Hamiltonian 
and Liouvillian dynamics. Using the Majorana basis de- 
fined above, one can write i? pare nt = c T Ei K ® l?) c = 
^c T Yc; the fluctuation matrix Y thus clearly appears as 
the counterpart of the parent Hamiltonian in the matrix 
representation given by Eq. (2). This one-to-one corre- 
spondence between H palCDt and Y clearly shows that the 
latter contains at least partial information about the dis- 
sipativc dynamics. As we now proceed to prove, the in- 
formation provided by H pareDt -or, equivalently, by Y- is 
sufficient to fully describe the dissipative dynamics when- 
ever the Lindblad operators satisfy {Li,Lj} = for all 
i,j. This result will emerge as a by-product of the more 
general result proved below, which can be stated as fol- 
lows: 

Theorem. For the purely dissipative dynamics defined 
by Eq. (1), the following statements are equivalent: 

(i) There exists initial conditions leading to a pure 
steady state; 



(ii) The dissipative dynamics admits an equivalent 
Hamiltonian description, and the latter is pro- 
vided by the parent Hamiltonian H palCDt = 

(iii) The Lindblad operators form a set of anticom- 
muting operators, i.e. {Li,Lj} = for all i,j 
€ {l,2,...,n}; 

(iv) [X,Y] = and X 2 = -\Y 2 . 

If there are as many non-zero Lindblad operators as 
fcrmionic sites -such that the Lindblad operators form 
a complete set of anticommuting operators- these condi- 
tions ensure that the steady state is pure and unique. 

Proof. Let us first assume that {Li, Lj} = is satisfied 
for all In the matrix representation defined above, 
this condition translates as iflj = = for all i,j. 
Remembering the explicit form of the X and Y matrices 
(see Eq. (3)), we then find 

[X, Y] = [X,iY] = and X 2 = -\Y 2 = (|F) 2 , (8) 

as can easily be verified. The first equation above im- 
plies that the Hermitian matrices X and iY can be 
diagonalized simultaneously by a unitary transforma- 
tion U. We denote the corresponding diagonal matri- 
ces as D x and D Y , and notice that D Y = iWYU = 
diag(Ay j i, — Ay,i, . . . , Xy,n, — Xy,n) with Xyj > must 
be satisfied owing to the antisymmetry of Y. Us- 
ing the second equation above, we then obtain D\ = 
\Dy, which in turn implies D x = U^XU — 
diag(Ax,i, A x ,i, • • • , X X . N , X x ,n) with X XJ = X Y ,j/2 > 
(remember that X is positive semidefinite) . The spec- 
tra of X and iY are thus in one-to-one correspondence 
(in particular, the spectrum of X is doubly degenerate). 
Therefore, the fluctuation matrix Y -or, equivalently, the 
parent Hamiltonian i? pa rcnt~ contains the same amount 
of information about the dissipative dynamics as the 
damping matrix X. In other words, the parent Hamil- 
tonian provides a complete description of the dissipativc 
dynamics when the Lindblad operators form a set of an- 
ticommuting operators. 

As dictated by Eq. (2), the steady-state correlation 
matrix T must satisfy the equation {X, F} = Y, which 
can equivalently be written as {D X ,T'} = D Y with f" = 
iU'TU. It should be clear that this steady-state equation 
does not determine the form of T' in the decoherence-free 
subspace corresponding to X X j = (as argued previ- 
ously, the latter solely depends on the initial conditions) . 
However, one can easily verify that the steady-state equa- 
tion constrains T' to a block-diagonal form (Bj \o v in 
the subspace corresponding to X x ,j > owing to the 
tight relation between D x and Dy derived above. With 
proper initial conditions, it is therefore possible to ob- 
tain f' = ejLi i°- y , which satisfies f' 2 = f 2 = -1. This 
concludes the first part of our proof. 
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Let us now show that the reverse is true, namely, that 
the existence of a pure steady state T implies {Li ,Lj} = 
for all We first notice that [f, Y] = [f,{X,f}} = 
directly follows from the fact that T satisfies {X, T} = Y 
and F 2 = — 1. ir and iY can thus be diagonalized simul- 
taneously by a unitary transformation U, yielding Df = 
iWTU = diag(l,-l,...,l,-l) and D Y = iWYU = 
diag(Ayi, — Ay,i, • • ■ , ^y,n, — Ay,jv) with Xy.j > 0. Ex- 
pressing the steady-state equation in the basis defined 
by U, we then obtain {fflXU, T'} = Y' , which in turn 
implies WXU = diag(A x ,i, Xx,i, ■ ■ ■ , Ax,jv, A x ,jv) with 
Xx,j = > 0. Similarly as above, the spectrum 

of X is therefore in one-to-one correspondence with the 
spectrum of Y; in particular, it is doubly degenerate. 
One can easily verify that X 2 = —\Y 2 , as in Eq. (8). 
Remembering the explicit form of the X and Y matrices 
(see Eq. (3)), we thus obtain 

£(i* ® if ® ij) = £(if i.-xi* ® ij) = o, (9) 

ij ij 

which, since L and Ij have support in different subspaces 
for i ^ j, can only be satisfied if {Li,Lj} = iflj = for 
all This concludes the second and last part of our 
proof. ■ 

The above discussion shows that the damping and fluc- 
tuation matrices X and Y are in one-to-one correspon- 
dence -and therefore both have full information about 
the dissipative dynamics- if and only if there exists a 
pure steady state or, equivalently, if and only if the Lind- 
blad operators form a set of anticommuting operators. 
Remembering that the fluctuation matrix is the exact 
counterpart of the parent Hamiltonian (i.e. i? pa ront = 
jC T Yc), we conclude that H parent fully characterizes the 
dissipative dynamics whenever the steady state of the 
latter is pure. In that case, its ground state exactly coin- 
cides with the steady state, and its ground-state degen- 
eracy -if any- reflects the existence of a decoherence-free 
subspace associated with zero-damping Majorana modes 
7 = 7^ satisfying the orthogonality condition {Li, 7} = 
for all i. 

In the more general case where the steady state of 
Eq. (1) is mixed, the one-to-one correspondence between 
X and Y breaks down and the parent Hamiltonian does 
not provide a complete description of the dissipative dy- 
namics anymore. The existence of Majorana zero modes 
is then solely determined by the form of the damping 
matrix X. In particular, a zero mode of H palCDt (or iY) 
does not necessarily correspond to a zero mode of X, 
although the converse is always true (see discussion fol- 
lowing Eqs. (4) and (5) above). Most crucially, any zero 
mode of -ff par cnt which does not coincide with a zero mode 
of X is effectively traced out of the system in steady state, 
thereby giving rise to an intrinsic purity zero mode. This 
effect lies at the very heart of the model constructed in 
the main text: in that model, H parent exhibits isolated 



pairs of Majorana zero modes in vortices of odd vortic- 
ity or on the edge of the system. However, one member 
of each pair is traded for a purity zero mode as a result 
of dissipation, so that single Majorana zero modes are 
obtained. 

To conclude, we remark that the well-known Hamil- 
tonian mechanisms of topological protection embedded 
in .ffparcnt do not extend to X. Indeed, while the anti- 
symmetry of Y ensures that the Majorana zero modes 
of -ffparcnt come in pairs and are thus unaffected by local 
perturbations when spatially separated, the symmetry of 
X does not allow for such a mechanism. Majorana zero 
modes arc therefore generally not protected against arbi- 
trary local dissipative perturbations (as shown explicitly 
above). Although this effect peculiar to dissipation is po- 
tentially harmful, as has been recognized in Refs. [5, 6] 
in another context, the model constructed in the main 
text shows that it can also lead to intriguing physics 
such as the transformation of Abelian vortices into non- 
Abelian ones. In general, engineered dissipation allows 
us to prepare topological phases with unexpected edge 
physics and no Hamiltonian counterpart. 

CALCULATION OF THE TOPOLOGICAL 
INVARIANTS AND RELATION BETWEEN 2D 
AND ID MODELS 

In the following section, we give an explicit derivation 
of the Chern and winding numbers presented in the main 
text (sec Eqs. (3) and (5)). We recall that both of these 
quantities are bulk topological invariants, which are only 
defined in the limit of an infinite system. In what follows, 
we will therefore always assume the system to be infinite 
and translation invariant, for convenience. To conclude 
the section, we provide the main ingredient for mapping 
the 2D vortex model to a ID model as argued in the main 
text. 

2D case: Chern number 1/2D 

Using the notations of the main text, we consider the 
general case of a 2D system driven by momentum-space 
Lindblad operators of the form Lk = Uk^k + V] l a < _ k , and 
define a "Fermi surface" F as the following set of points 
in the Brillouin zone (BZ): 

^={k: |^k| = l}, (10) 

where tpk = Icicle 1015 is the so-called pair wavefunction, 
defined as </?k = v k/ u k- Since the Brillouin zone cor- 
responds to a torus, this Fermi surface does not nec- 
essarily admit a well-defined interior and a natural ori- 
entation. Following Rcf. [7], we define an electron-like 
region £ = {k : |</?k| > 1} and a hole-like region 
W = {k : |<pk| < 1}, and choose the Fermi surface to 
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be oriented positively whenever the electron-like region 
is on the left as we go along it. With this convention, 
the interior of T can be identified with £. In general, 
the electron-like region corresponds to disconnected (sim- 
ply or multiply connected) regions of the Brillouin zone, 
and the Fermi surface JF can be described as a finite 
union of piecewise smooth, simple closed curves F\, i.e. 
T = U A J~\ (more pathological Fermi surfaces can be en- 
countered, but the Chern number can only be defined 
in that case). As shown in Rcf. [7], the Chern number 
can then be expressed as a sum of winding numbers W\ 
around each of the closed curves F\, namely, 



^2D = 7- / cPk n k • (9 fcx n k x d k n k ) 
4tt J bz 



(ii) 



with winding numbers 



2tt 



V k 6» k • dk 



= h £ + d ky kdk y ). (12) 

If du dk„ # k and dk„ dk, r # k arc continuous functions on the 

j; y y 

electron- like region enclosed by F\, Green's theorem can 
be applied and gives 

— f f (dk^d^Ok - d ky d k J k ) dk x dk y = 0. 

Z7r J JintFx 

(13) 



Wx = 



In other words, W\ trivially vanishes. In order to allow 
for non-trivial contributions to the Chern number, it is 
therefore necessary to have discontinuities in the second- 
order mixed derivatives of # k in the electron-like region 
(as a matter of fact, in the absence of such discontinuities, 
one can always perform a gauge transformation to get rid 
of the phase # k ). We note that a necessary and sufficient 
condition for a particular discontinuity point kj to be 
contained in £ is given by 



lim 



(|¥>k|-l)>0. 



(14) 



When this condition is satisfied, kj contributes to the 
Chern number by an integer value Wi equal to the num- 
ber of times that the phase # k winds by 2ir when going 
around kj (in the counter-clockwise direction), as dic- 
tated by Eq. (12). This provides us with a simple way of 
determining the Chern number. 

Let us now examine the specific case considered in the 
main text, with Lindblad operators defined as 



Li = pa\ + (a* + a* + ot_. + at.) 



+ (a, 



+e x 



lfl 



i+e y 



ia»-e y ), 



(15) 



where we have set a = 1 and </> = (see main text), with- 
out loss of generality. Expressed in momentum space, 
these operators take the form 



-«ka^ k , 



l-k = "k»k , - 1, 

u k = 2i(sin (k x ) + isin (k y )), 
"k = 8 + 2(cos (k x ) + cos (k y )). 



The minimum of their "norm" 



7V k = {4,i k } = |w k | 



(16) 



(17) 



can easily be seen to be proportional to the dissipative 
gap. Since it vanishes for — — 4 at (k x ,k y ) = (0,0), 
= at (k x ,k y ) = (tt,0), (0,tt), and = 4 at (k x ,k y ) = 
(it, it), we anticipate the existence of four topological 
phases, depending on the value of 0. As shown above, 
the Chern number is solely determined by the existence 
of discontinuity points and the behavior of the pair wave- 
function ipk around them. Here we have 



</?k = 



+ 2(cos (k x ) + cos (k y )) 
2i(sin (k x ) + i sin (k y )) 



(18) 



and the phase k is defined everywhere except at the 
four discontinuity points kj € {(0, 0), (0, tt), (tt, 0), (tt, tt)} 
which coincide with the gap-closing points found above 
(note that these points are invariant under time-reversal 
symmetry kj — > — kj). The behavior of <y? k around them 
is given by 

+ 2 [cos (e cos 6) + cos (e sin 0)1 

lim (£> k = lim : — 

k->(o,o) c^o 2i [sin (e cos 6) + isin (e sin 9)\ 

= lim e 

e^o 2ie 

+ 2 [cos (e cos 9) + cos (tt + e sin 0)1 

lim (/? k = lim - — 

k^(o,7r) e^o 2i [sin (e cos 9) + i sin (tt + e sin 9)\ 

= lim ^e w 
e^o 2ie 

+ 2 [cos (tt + e cos 9) + cos (e sin 9)] 

lim tp k = lim — . — -. — — — - — : — — T ± 

k^(7r,o) e^o 2i [sin (tt + ecosfl) + ism (esmft)J 

e-^*6 -2ie^ 

+ 2 [cos (tt + e cos 9) + cos (tt + e sin ( 

lim i/? k = lim 

k^(7r,7r) c^o 2i sin (tt + e cos 6) + i sin (tt + e sin 6) 



= lim 



8-4 



e-s-o -2ie 



(19) 



where e, 9 are polar coordinates defined around the point 
of interest. The phase 9k is thus found to rotate in 
the counter-clockwise direction around the discontinuity 
points ki = (0,0) and k 2 = (tt, tt), and in the clockwise 
direction around k 3 = (0,7r) and k4 = (tt,0). In light 
of the above discussion, ki and k 2 (k 3 and 1*4) therefore 
both contribute to the Chern number by +1 (respectively 
— 1) if Eq. (14) is satisfied. As can easily be seen from 
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Eq. (19) above, we thus obtain i^d = +1 + 1 — 1 — 1 = 
for all values of (3 except at the isolated points (3 = and 
(3 = ±4 where the dissipative gap closes. 

We finally indicate the symmetry class of the 2D bulk 
system according to the Hamiltonian classification of 
Ref . [8] . Such a classification is made possible here due to 
the fact that the steady state of the translation-invariant 
system is pure, and thus identical to the ground state of 
the parent Hamiltonian H parcnt . Since the latter only 
has particle-hole symmetry (C) with unitary operator 
Uc — X , the 2D model belongs to the D class. 



ID case: winding number 

We now examine the case of an infinite, translation- 
invariant ID system with Lindblad operators as defined 
by Eq. (4) in the main text, namely, 

U =pa\ + (at.! + ot +1 ) + (-Oi_i + Oi+i), (20) 



or, in momentum space, 



L k = 2i sin {k)ak + (/?' + 2 cos {k))a 



(21) 



As in the 2D case discussed above, the dissipative gap 
is proportional to the minimum of the norm A4 = 
{Li, L k } = j3' 2 + A(3' cos (fc) + 4. Since the latter van- 
ishes for j3' = ±2 (at k = ir and k = 0, respectively), 
we must distinguish three regions in parameter space, 
namely: (3' < -2, |/?'| < 2, and (3' > 2. As long as 
| 7^ 2, we can define normalized Lindblad operators 
L k L k /^/J\fk = u k cLk + Vk<r_ k , with 



u k = 



2i sin (k) 



Vk 



^/3' 2 +4(3' cos(fc) +4' 

13' + 2 cos (fc) 
a//3 ,2 +4/3' cos(fc)+4' 



(22) 
(23) 



and represent the latter by a real unit vector n k = 
( n fc' n fc' n fe) = (0, — 2m^ Imw^,, u 2 — l u fc| 2 ), which can be 
cast into the form of a complex number e l ^ k = n k + \n v k 
for some angle <j> k . Based on this construction, the so- 
called winding number topological invariant then takes 
the form [9, 10] 

v m = f dk&- (n fe x d k n k ) 

27T JBZ 



1 f J, / z dn l 

— / dk(n z k —f- 
271" Jbz dk 



k dk ' 



l 



7T 



dk 







<PlT - 00 



7T 



d(j>k 
Ilk ' 

e z. 



l 



7T 



(24) 



The associated topological phase diagram can be derived 
by examining the behavior around specific points located 



in the three parameter regions (3' < —2, \(3'\ < 2 and (3' > 
2 where the calculation is most simple and the topological 
features easiest to visualize. Here wc find 



13' -> -co : u k = -1; v k = 0; 

n k = (0,0,1); cp k =0, 
= : Mfc = cos (fc); Vfc = — isin(fc); 

rife = (0, sin (2k), cos (2fc)); 0^ = 2k, 
13' -> +oo : u k = 1; v k = 0; 

n fc = (0,0,1); 4> k =0. 



(25) 



As given by Eq. (24), we thus obtain u\y> — 2 for \f3'\ < 2 
and vib = otherwise. 

We finally remark that the ID model is chiral and be- 
longs to the BDI class. Using the fact that the bulk 
steady state is pure -as in the 2D case above- this can 
easily be shown by examining the symmetries of the asso- 
ciated parent Hamiltonian. The latter has time-reversal 
symmetry (T) with unitary operator Ut = 1 (which triv- 
ially satisfies U^Ut = +1) and particle- hole symmetry 
(C) with unitary Uc = & x (which satisfies U c Uc = +1), 
and thus indeed belongs to the BDI class. 



Mapping the 2D vortex model to a ID wire 

We now provide further details regarding the mapping 
to a chiral ID wire presented in the main text. For sim- 
plicity, we work in the continuum limit and set a = 1, 
<j) = in the general form of the cross Lindblad opera- 
tors (Eq. (2) in the main text). In an annular region 
centered around the vortex core (with r > r c ; see main 
text), these operators take the form L(x, y) = ((f3 + 4) + 
dl + d 2 )tp^(x,y) + e-^(d x + id y )ip(x,y) where ip(x,y) 
denotes the fermionic field and (x, y) are Cartesian co- 
ordinates on the plane. Expressed in polar coordinates 
(r, if) defined from the center of the vortex core, they be- 
comc L(r,<p) = {{(3 + 4) + 0? +d*/r 2 + (l/r)d r ) V> f (r, ip) + 
(d r + \d v /r) ip(r, <p), and the external U(l) gauge field 
e~ %Lp describing the vortex phase does not appear explic- 
itly anymore. The annular region can be mapped onto 
a cylinder with Cartesian coordinates (x',y') by defin- 
ing x' = r, y' = rip (see main text, Fig. 2(a), right). 
The Lindblad operators then transform as L(x',y') = 
(08 + 4) + d 2 xl + d 2 ) ^{x', y>) + [d x ,+ id y ,)ip{x', y') [11], 
and we recover the form of Eq. (2) in the main text cor- 
responding to our original translation-invariant model on 
the plane. This shows that there exists a one-to-one cor- 
respondence between our original model embedded in the 
plane with a single £ = 1 vortex and the same model in 
cylinder geometry with no vortex (i.e. with p-wave oper- 
ator d x i + idyi). Physically, this stems from the fact that 
the gauge field e~ tv that was imposed on the plane to de- 
scribe the vortex naturally arises on the cylinder owing 
to the extrinsic curvature of the latter. 
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Starting from the cylinder model, one can perform a 
Fourier transform in the translation-invariant direction 
of the cylinder in order to decouple the 2D model into 
a stack of ID wires. Two of them, associated with the 
modes corresponding to the momenta and ir, turn out 
to exhibit parameter regimes in which they are topo- 
logically nontrivial, as indicated by the winding number 
topological invariant (see main text). 



MAJORANA ZERO MODES IN THE CHIRAL ID 
WIRE 

We consider the extended ID wire depicted in Fig. 1(c) 
in the main text, with position-dependent Lindblad op- 
erators of the form 

Li = ft at + aj_ 1 + a\ +1 - /(ri_i)aj_i + f(r i+1 )a i+1 , 

(26) 

where f(r) denotes the vortex profile function (see main 
text) and r the position along the wire (this notation em- 
phasizes the fact that the ID wire effectively describes the 
behavior of a dissipative vortex with vorticity £ = 1 in the 
radial direction away from its core, as discussed in the 
main text). Since we are mainly interested in the behav- 
ior of the original 2D system inside the vortex core, we 
assume that the wire is semi-infinite, with open bound- 
ary conditions on its left side imposed by the geometry 
of the vortex core alone, i.e. by the profile function f(r). 
In the framework developed in the first section above, 
the Lindblad operators take the form of complex vectors 
li (defined such that L, = if c where c is a vector of 
local Majorana operators; see above) with non-zero com- 
ponents 

(Ii) 2 (i-i)-i = (1 - /(ri-i))/2 = <5(ri_i), 
(li) 2 (i-i) = i(-l - /(ri-i))/2 = -i(l - <5(ri_i)), 
(l*)2i-l = 072, 
(khi = -i/372, 

(I*) 2 (i+i)-i = (1 + /(r i+ i))/2 = 1 - S(r i+1 ), 
(Ii) 2 (i+i) = i(-l + /(r i+ i))/2 = -i*(r i+ i), 



(27) 



where S(ri) is the deviation from the value /(r») = 1 
found in the bulk, defined such that /(r^) = 1 — 2<5(rj) 
with < S(ri) < |. As shown previously, a real unit 
vector v defines a Majorana zero mode 7 = v T c if and 
only if if v = for all i. Noticing that the components 
are here real (respectively purely imaginary) for odd 
(even) indices j, we can restrict our search for Majorana 
zero modes to vectors v of the form 

(i) (v) 2j _ 1 =a j ,(v) 2j =0, (28) 
or(ii) (v) 2j _! = 0, (v) 2j = aj. (29) 

with aj e BL As will be argued later, both cases give 
similar results owing to the symmetry of Eq. (27) when 



considering even and odd components separately. We 
thus focus on (i), in which case the necessary and suffi- 
cient condition if v = translates as 

a j+ i(l - <5(r<+i)) + on{0 12) + a<_i<f(ri_i) = 0. (30) 

In general, such a homogeneous linear second-order re- 
currence relation cannot be solved in closed form since its 
coefficients are not constant. In order to make progress, 
we examine the possible solutions away from the vortex 
core (i.e. away from the left edge of the ID wire) where 
the vortex profile function is constant, namely, 6{r{) = 5. 
In that case, one easily finds the general solution 



a. 



c + (r + y + c_(r_y 



with C±€l and 



-P'/2±y/(/3'/2)*-4S(l-6) 
2(1 -S) 



(31) 



(32) 



Since r± must be real, this solution exists if and only if 
I > 4-^/(5(1 — 5), which is trivially satisfied away from 
the vortex core where /(t») = 1, i.e. 5 — 0. The general 
solution then reads a, = C(— fi' /2) % with C E K, and 
is normalizable if and only if \/3'\ < 2. Assuming that 
this condition is satisfied, the solution can be extended 
inside the vortex core using the recurrence relation of 
Eq. (30), provided that the extension is consistent with 
the boundary conditions imposed at the center of the 
vortex core (i.e. on the left edge of the wire). Owing 
to the rotational symmetry of the vortex, the boundary 
conditions must be defined as 



(33) 



where we have chosen the index of the leftmost site of 
the wire as i — 0. Using Eq. (30), we then find ct\ — 
Oio{—0'/2). We must therefore have ao > in order to 
obtain a non-trivial solution. Since v must be normalized 
to unity, there is no freedom in choosing the value of a a , 
and the solution of Eq. (30) is unique. This proves that 
there exists a single Majorana zero mode localized in 
the vortex core (or, cquivalcntly, on the edge of the ID 
wire) in the whole parameter range < 2, and that the 
latter decays exponentially away from the vortex core, on 
a characteristic length scale £ = - l/|log (|/3'|/2)| 

which diverges as |/?'| — > 2, as argued in the main text 
(and as verified numerically). 

Remembering that the Chcrn number vanishes for 
\(3'\ < 2 (i.e. for < \f3\ < 4; see main text), the ex- 
istence of a single Majorana zero mode in the vortex 
core is somewhat surprising: according to the well-known 
bulk-edge correspondence for the Chcrn number [8, 12- 
14], one would naively expect the existence of an even 
number of such modes. We wish to emphasize that 
there is no contradiction here, since bulk-edge correspon- 
dence arguments only apply to the parent Hamiltonian 



8 



-ffparent = £\ ijij (or to the matrix iK) defined above; 
not to the damping matrix X which solely determines the 
existence of Majorana zero modes in our dissipative set- 
ting (see Eq. (3) and discussion thereof). In the present 
case, the bulk spectrum and eigenmodes of -ff pa rent ex- 
actly coincide with those of X due to the fact that the 
Lindblad operators anticommute with each other in the 
bulk. However, the fact that these anticommutation rela- 
tions are not satisfied in the vortex core leads to a loss of 
purity. In light of the general discussion presented in the 
first section above, we expect the bulk-edge correspon- 
dence pertaining to H paient = J2i ^l^i to be reflected in 
the appearance of a purity zero mode, as seen numerically 
(see main text). 

We finally remark on the case given by Eq. (29) which 
has not been considered so far. Noticing in Eq. (27) that 
the odd components (lj)2j-i ~for increasing j- have the 
same form as the even components (l,)2j ^for decreasing 
j- one can easily convince oneself that Eq. (29) corre- 
sponds to the symmetric situation of a single Majorana 
zero mode localized on the right side of the ID wire. 
This can be shown explicitly by imposing open dissi- 
pative boundary conditions on the right of the ID wire 
(thereby introducing an edge in the original 2D system) 
and repeating the construction above. Since such dissi- 
pative boundary condition appearing on (and defining) 
the physical edge of the system far away from the vor- 
tex core are mathematically similar to those appearing 
in the vortex core, a single Majorana zero mode must 
also be found on the edge of the system. In general, 
these boundary conditions will not be compatible with 
the Majorana zero mode localized in the vortex core; 
similarly, the boundary conditions imposed in the vor- 
tex core will not be compatible with the Majorana zero 
mode found on the edge of the system. As a result, both 
Majorana zero modes will acquire a small damping rate 
A ~ (\(3'\/2) 2N which decreases exponentially with the 
system size N, i.e. with the distance between them. The 
situation is therefore reminiscent of that of interacting 
Majorana modes in a Hamiltonian system: \j3'\ can be 
viewed as a coupling constant and A as a small splitting 
of the Majorana modes caused by inter-edge tunneling. 

IMPLEMENTATION IN ULTRACOLD ATOM 
SYSTEMS 

Here we discuss how to implement the dissipative dy- 
namics investigated in the main text (see e.g. Eq. (2)) 
with spinless fermionic atoms in optical lattices. In ultra- 
cold atom experiments, spin is not a fundamental quan- 
tity, but is realized by a specific hypcrfine state. If only a 
single hyperfine level of the atoms is populated, one ob- 
tains an ensemble of fermions that are effectively spinless. 
This underlies the possibility of realizing the Kitaev's 
Majorana wire with cold atoms [15]. The idea of the 



implementation of the situation with no vortex follows 
Rcfs. [10, 16]; we summarize its main ingredients below. 
We first sketch how microscopically number-conserving 
Lindblad operators, corresponding to a quartic Liouvillc 
operator, relate to the quadratic master equation of the 
main text in the long-time limit. We then show how to 
implement the number-conserving operators in practice, 
and in particular how to optically imprint a dissipative 
vortex as defined in the main text. 



From number-conserving quartic to quadratic 
dissipative dynamics 

We start from a dissipative dynamics that microscop- 
ically conserves particle number, 

N 

dtP = ko Yl { e iP £ l - P}) > ( 34 ) 

i=l 

with bilinear Lindblad operators of the form £j = CjAi, 
where Ai (Cj ) is a quasi- local superposition of fermionic 
annihilation (creation) operators. As anticipated in the 
main text, linear Lindblad operators such as the ones 
of Eq. (2) (see main text) then arise in the long-time 
and thermodynamic limit. More precisely, within a dis- 
sipative analog of BCS mean-field theory controlled by 
the proximity to the exactly known steady state (see the 
appendix of Ref. [10] for a detailed discussion), one can 
show that 

ti = CjAi ^ Li =Cj + ae^A t (35) 

for spinless fermions. This results in the following mean- 
field master equation governing the late-time damping of 
excitations close to the steady state: 

N 

d tP = K Y, { L *P L l - U L \Li,p}) ■ (36) 

i=l 

In Eq. (35), the phase <j> originates from spontaneous U{1) 
symmetry breaking, i.e., it is not fixed by the dissipative 
dynamics [17]; it will be set to zero in the following, with- 
out loss of generality. The relative strength a > be- 
tween creation and annihilation parts, on the other hand, 
is fixed by the system filling according to a specific num- 
ber equation (see Eq. (39) below). For the homogeneous 
setting in the thermodynamic limit, this number equa- 
tion -as well as further characteristics of the late-time 
dissipative dynamics- can be specified explicitly. Using 
similar notations as before (see e.g. Eq. (16)), the Fourier 
transform of Eq. (36) can be written as 

dtP = J (0^ k [L kP Lt p}) , (37) 
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with [18] 



£k = W k 1/2 L k , L k = ("kak + awka^ k ), 



A4 



|a«k| 2 , «k = K-A/k- 



(38) 



The explicit form of Uk, «k relevant to our model is given 
by Eq. (16). Note that we have kept the coefficient a 
explicitly (as opposed to setting a — 1 as in Eq. (15) for 
example) in order to examine how it relates to density. 
The fact that the pair wavefunction tp^ = f k /itk is anti- 
symmetric under momentum reflection (i.e. </? k = — ^>- k ) 
implies that {Zk, Z q } = {L^, L^} = for any pair of 
momenta k, q. In addition, the normalization chosen 
in Eq. (38) guarantees that {Zk,L q } = <5(q — k). The 
normalized Lindblad operators are thus proper fcrmionic 
quasiparticle operators; they generate a full Dirac algebra 
for the system and determine the steady state uniquely. 

In the homogeneous setting, the mean-field calculation 
yields a simple form for the number equation, 



n : 



/ 



d 2 q 



+ \av 



(39) 



which explicitly relates the parameter a to a given fixed 
average particle number. The effective damping rate en- 
tering Eq. (36) can also be calculated explicitly; it relates 
to the microscopic one of Eq. (34) as 



K = Kq 



BZ 



d 2 q \UqVg\' 

(2tt)2 | Uq 



av a 



(40) 



The simplifications arising in the long-time limit of 
Eq. (34) are similar to the ones emerging in the low- 
energy limit of a Hamiltonian theory describing a su- 
pcrfluid. Indeed, the low-energy physics of a fcrmion 
theory with weak attractive interactions (i.e. a quartic 
fcrmion theory) is universally described by a quadratic 
BCS Hamiltonian. In our dissipativc setting, the analog 
of the gap for single-particle excitations present in such 
theories is a dissipative gap R > 0, defined as the min- 
imum of the damping function K k , i.e. R = min k {Kk}- 
In the specific model considered in the main text (sec 
Eq. (2)), the damping function features a gap for all 
values of ft except for ft = 0, ±4. The dissipative gap 
closes at k* = (0, it), (tt,0) for ft — 0, k* — (ir,n) for 
ft = 4, and k* = (0,0) for ft = —4; in the vicinity of 
these points, the damping rate vanishes quadratically, 
i.e. K k ~ (k — k*) 2 . Physically, the existence of a dissi- 
pative gap (i) implies the exponentially fast convergence 
of all many-body observables to the steady state, and 
(ii) gives rise to a separation of time scales between the 
rapidly damped modes of the bulk and the (Majorana) 
zero-damping modes which are present for appropriate 
boundary conditions. 



Implementation of the number-conserving bilinear 
Lindblad operators 



The operators i t = CjAi appearing in Eq. (34) can be 
engineered via two-step processes involving an auxiliary 
degree of freedom, as shown in Refs. [10, 16, 19]. Auxil- 
iary degrees of freedom can be provided by the intermedi- 
ate sites of an optical superlattice, giving rise to a second, 
excited Bloch band (the target system residing in the 
lower sites of the lattice corresponds to the lowest Bloch 
band). The first step consists in coherently coupling a 
quasi-local superposition of fcrmions to the auxiliary de- 
gree of freedom, and is associated to the annihilation part 
Ai. The second step, on the other hand, introduces dissi- 
pation in the form of spontaneous phonon emission, and 
relates to the creation part C\: a decay channel for the 
particles from the excited Bloch band back into the low- 
est one is opened by coupling to a reservoir into which the 
whole lattice system is immersed. A suitable reservoir is 
provided by a bosonic supcrfluid; the bosonic nature of 
the environment is particularly useful, as we will see be- 
low. Microscopically, the system-reservoir coupling can 
be obtained from standard s-wave density-density inter- 
actions ~ fisys(a;)^bath(^) between system and bath par- 
ticles. Note that the resulting scattering processes con- 
serve the number of particles in the system. Nevertheless, 
the release of system energy necessary to decay sponta- 
neously from the excited to the lowest Bloch band is only 
possible via the emission of phonons into the superfiuid 
reservoir: energy exchange with the reservoir therefore 
occurs, while particle exchange does not. Ultimately, the 
coherent coupling to the auxiliary sites can be operated 
in such a way (far detuned) that the latter can be inte- 
grated out. We then end up with an effective dissipative 
dynamics for the lowest Bloch band only, described by 
operators ti = C[Ai -the annihilation part Ai results 
from the coherent coupling to the auxiliary site, and the 
creation part C\ from the spontaneous decay back to 
the target Bloch band. For a detailed exposition of this 
physics, we refer to [20]. 

We now examine how to implement the specific anni- 
hilation and creation parts relevant to our problem. Ai 
is obtained from coherent driving, i.e. Ai — J2j^ij a j' 
where Qjj is the Rabi frequency connecting to the aux- 
iliary site. This quantity is thus fully determined by the 
properties of the coherent drive laser. In particular, the 
short-distance p-wave structure of Ai can be realized with 
specific wavelength ratios and relative phases between the 
laser generating the optical lattice and the drive laser. 
The longer-range vortex profile, on the other hand, can 
be implemented with a drive laser carrying an optical vor- 
tex; such vortices are easily obtained e.g. in the far-field 
regime of a fork hologram or a spiral phase plate [21-23]. 
This results in an annihilation part with Rabi frequencies 



10 



of the form 

^ij ~ f( r ]) ell!V:i ( S i-e x ,j - &i+e x ,j + K S i-e y ,j ~ ^+e yJ )) , 

where f(r) is a slowly varying profile function propor- 
tional to the laser amplitude encoding the optical vortex, 
e.g. for a Laguerre-Gauss beam f(r) <~ {r/l)e~( r /^ with 
~> a (a being the lattice spacing). In this way, the 
optical angular momentum of light is directly imprinted 
onto the matter system. 

In contrast, the creation part results from spontaneous 
emission from the auxiliary site, and is therefore generi- 
cally isotropic. This gives C\ the desired s-wave symme- 
try (see Eq. (2) in the main text). The spatial extent of 
C\ , as well as the relative strength of the phase-transition 
parameter /?, can be tuned through the bath proper- 
ties [16]. With Ai = Qijdj encoding the vortex phase 
and profile, the form of the Lindblad operators associated 
with a dissipative vortex therefore precisely matches the 
one introduced in the corresponding section on p. 2 in the 
main text; it reduces to the desired translation- invariant 
form of Eq. (2) (main text) for constant Rabi frequency. 
We have verified numerically that the precise structure 
of the vortex (different choices of f(r)) does not affect 
the findings presented in the main text. 

Finally it is clear that configurations of several vortices 
can be produced optically, and that single vortices can in 
principle be moved adiabatically around the lattice by 
displacing the optical device generating the optical vor- 
tex. Since (i) adiabatic parameter changes of the Liou- 
ville operator give rise to a gauge transformation just as 
in a Hamiltonian dynamics [10] and (ii) the vortices typ- 
ically carry single (unpaired) Majorana modes (see main 
text), the exchange statistics must be non-Abelian [24]. 

Discussion 

We comment on the conservation of particle number 
and parity in our system. A first key point of our set- 
ting is that it relies microscopically on the dissipative 
coupling to a bosonic environment -in contrast to many 
solid-state implementations of fcrmionic superfluids or 
superconductors- such that the fermion parity of the tar- 
get system is an excellent quantum number: any coupling 
of bosonic degrees of freedom to fermions must be even 
in the number of fermion operators due to superselection 
rules. This is an important prerequisite for robust topo- 
logical order [25, 26]. Second, the coupling to the reser- 
voir is such that the fermionic particle number is exactly 
conserved microscopically (that is, [£i,N] — for all i, 
where N = . is the total fermionic particle num- 
ber), which obviously implies exact parity conservation. 
Exact number conservation is given up subsequently in 
the mean-field theory, in full analogy to the BCS strat- 
egy and justified by the equivalence of fixed-number 



and fixed-phase BCS states in the thermodynamic limit, 
where relative number fluctuations in the fixed-phase 
BCS wavefunction scale as <~ 1 /VN (N being the number 
of degrees of freedom). While it may seem as if particles 
are lost or created in the system due to the interaction 
with the environment from the linear Lindblad opera- 
tors, this must really be understood as processes that 
annihilate or create fermions into or out of the fermionic 
pair mean field. This is particularly transparent from the 
explicit form of the effective mean-field damping rate k 
of Eq. (40), which includes weighted integrals of corre- 
lation functions such as (a q a_ q ) <~ UqVq providing for 
a mean field for the mode k (see Refs. [10, 17] for a de- 
tailed discussion). In particular, this entails that the par- 
ity of the fermionic target system is fixed via its initial 
state. Finally, we critically examine the boundary con- 
ditions in this light. The vortex is generated optically, 
and crucially does not modify the property of particle 
number conservation of the Lindblad operators. There- 
fore, no problem arises in vortices with particle or parity 
conservation. Particle number conservation at the edge 
of the system, in contrast, is a subtler issue. It could 
be achieved along the lines of Ref. [15], but the need for 
such subtleties can be avoided by trapping any Majorana 
mode localized on the edge in an additional vortex. 
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